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INTRODUCTION 

The present investigation was carried out to determine the efficiency 
and accuracy of a moment method as applied to the calculation of incompres- 
sible turbulent boundary layers. This method, at least conceptually, promises 
to be more accurate than the classical von Kazaas integral methods while 
requiring leee computing tine than previously employed finite difference 
methods. 

The moment method employed In the preseat study is popularly known 
as the ’’method of Integral relations" , following Dorodnitsyn ( 1960 ) . It 
is one of a family of methods that are kaouo in the sathexatics literature 
aa Petrov's generalization of the Galervdn r ethod (Xlkhlin and Saolitskiy 
(1967)). Briefly, the mathematical procedure employed is to approximate 
the actual solution to a differential equation, where the solution is given 
by aa infinite series (hence, the solution i said to lie in an infinite 
dimensional vector space). The appro* imatit a is taken in the fora of a 
partial sum of the Infinite series, the 3 placing the approximate solution 
in a finite dimensional vector space. Vais finite dimensional space is 
completely described by a system of coordinate vectors (normally referred 
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to in moment methods as "weighting functions") in that any decent (e.g. 
th« assumed solution) nay be expressed as a linear combination of the coordi- 
nate vectors. The dimension of the subspace chosen to represent the solu- 
tion is sons tines known as the "order of approbation", and is equal to the 
aueber of paraMters which nay be determined from use of the moment method. 

The various nome n t methods used for solving the boundary layer equations 
basically differ only in the chitci of these coordinate vectors and the 
form ass umed for the approximate solution. The success of any moment method, 
then* depends on how well these choices can represent the actual solution. 
Several constructive theorems demonstrating exactness and convergence of 
general moment methods applied to non-linear equations have been given by 

s 

Fetryshyn (1967). 

Moment methods similar to that employed herein have been employed in 
the calculation of incompressible laminar boundary layers by Dorodnitsyn 
(1960) and Bethel and Abbott (1966) and to compressible laminar boundary 
layara by Pallooe (1961) * Pavlovs ltii (lSf2) , Nielsen et al (1965) , and 
Holt (1966). The adequacy of moment methods as applied to laminar flows 
has been demo ns trated by these studies. 

Deiwert sad Abbott (1967) used a moment method with a three parameter 
velocity gradient profile together with empirical shear integrals, representing 
the turbulent transport mechanism, to solve the incompressible turbulent 
boundary layer equations. Through the use of empirical shear integrals 
determined from selected sets of data, they were able to predict the flow 
quantities for those particular sets of data* thus demonstrating the adequacy 
of the moment method applied in turbulent boundary layer flows. 
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In the present investigation, the results obtained using two, three, 
and four parameter velocity profiles were compared with flat plate data. 

It was found that the accuracy improved with Increasing order of approxima- 
tion. Based on these flat plate studies, it vas felt that a four parameter 
velocity profile provided sufficient accuracy, and that higher parameter 
profiles would require excessive computing time. All of the results obtained 
la this study employed the four parameter profile. In order to ;rovlde a 
general procedure for predicting the behavior of the turbulent boundary 
layer with arbitrary free stream velocity variation, an eddy viscosity 
model vas employed for the numerical evaluation of the shear integrals. A 
computing program was generated in this investigation and has been written 
so that any model of the turbulent transport mechanism may be easily incorpor- 
ated. The results of this program are compared with the experimental data 
uf sixteen turbulent flows. These flows are listed in table 1. 



NKiEXCLATU'Rl* 


Cj(s) coefficients In approximate velocity gradient profile 

representation 

f^(0) »g^(0) weighting functions employed in equations (5) and (6) 
INT 

i dissipation Integrals defined in equation (10) 

P,Q,R exponents defined in equation (8) 

s transformed 6treaavise coordinate 

u C/U. 

a exponent defined in equation (10) 

X eigenvalue of linearized version of equation (10) 

a 2 error magnitudes defined in equations (3a) and (4a) 

0 reciprocal normalized velocity gradient, «y/3u 

Superscript# 

* denotes approximate functional representation 

Subscript 

o denotes quantities at initial station 


*The uniform symbols suggested by Kline et al have been employed wherever 
possible. This table defines only those symbols which are not provided 
for in the uniform symbol tabulation. 
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The analysis given below defers in concept and procedure from those 
previously published. It is the authors' belief that the present analysis, 
vhilt more complex mathematically , provides a clearer understanding of 
the operations employed, than do previous derivations. 

The boundary layer equations for two-dimensional incompressible turbulent 
flow are written in Boussinesq form as: 


dy 


XJ m 




( 1 ) 




( 2 ) 


where •• <5^ 


The terms and have been considered negligible in comparison 

with c>eix > . 

ay 

Employing a change of variable: u ■ U/U^, 

the Crocco transformation x— *-s(x), y — *-u(x,y), 

and setting 9 * ^ 

®y 

ve obtain from (1) and (2) : 
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tern; ~r~ is 

o S 


given by 


^ f ft cr M. 


from the definition of 0/. 


As this point we have two partial differential equations in the two unknowns 
V and 0. These equations differ fror. (1) and (2) in that now one of the 
independent variables, (i.e. u) assumes values in the finite interval 0,1 

We now employ Kantorovich's technique to reduce the partial differential 
equations to ordinary differential equations by approximating the u dependent - 
of the unknowns V and 0. This approximation is accomplished by the generalised 
Galerkin procedure (see Kantorovich and Krylov (1964) and Kikhlin and 
Scolitskiy (1967)). In this procedure, it is assumed that the functional 
dependence on u of the unknowns V and 0 can be represented by partial suss 
of the form: n - n * 

■ £ as % a ^ C I C3> u. 1 " 


respectively. 

The significance of the transforms employed to obtain equations (3) 
and (4) now becomes evident. From the Teiers trass approximation theorem, 
it is clear that the chosen form of V and y can be made to respectively 
represent V and 0 as accurately as one desires on the closed interval 0,1 
If we now substitute the approximate values of 0 and V, i.e., 0 and V, 

into equations (3) and (4) we obtain: 

•» 

IVk 

?du + v ^ Z Ci- O'idlir*, _x A r(\ + *>\-v 1 * <r K 


(3a) 





ffu ctu^ 4 , J. 


(4a) 
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where o^ and differ from /.er-i by virtue of the fact that approximate 
representations of 0 and V have been introduced. We now require the errors 
and o^ to be orthogonal to two linearly independent sets of coordinate 
vectors denoted by "f and respectively. The form of and 

are given by: 

%*• _ y _ 

where B and D are arbitrary linear operators (c.f. Kikhlin and Snolitskiy 
(1967)). Invoking the infinite dimensional form of the inner product to 
Impose the orthogonality condition over the domain 0 <_ u 1, ue obtain: 

A U ** u 

\ ^ da + § - £ (5 > 
L u « ^ do> g r ] ^ 

O »0 w ^ J 




1 Vv. f- f j_ $ u dVc p 4 . _L cS v 1 clu - o 

^0*0^ ^ Uco ^ ^ 


VO) 


<d S t) 6 

The inner product employed is a bi-linear operator so that (5) and (6) 
may be added and, upon carrying out the indicated operations and requiring 
that g. - ^Ji w and f. (1) * 0, written as; 

. 1 ^ . \ du*r r^o + -\§-vi 

* r rr 31 \ ad'* uJ ?-yd /) 

J<v, A *;U 


U *\ 


•S.J. iueswi d “ - 

**u -c 


where primes denote differentiation with respect to u. 




The choice of ^(u) is noc arbitrary, but, rather, is restricted to 
functions satisfying the boundary conditions @ u “ 0 and u • 1. One such 0 is 
that which has been employed in laminar boundary layer studies: 

3 *_L- | . 

* <-D-£. 4 <*> 

The form of the assumed velocity gradient profile, employed in the present 
study is : 

^ i-Ct 

where P,Q, and R are exponents assumed to be aibitrary for the moment. 

. The fj(u) ; chosen in the present study are: ' 

1 i - 1,2, 3, 4. (9) 

Substitution of (8a) and (9) into equation (7) and integration of the u 
dependence yields a system of 4 ordinary non-linear differential equations 
in the four unknowns Cj(s). These equations may be written in matrix form as: 


Wi XKT ' 


where 

a 

for 


'Aco^l 


.. {P* L{1 

'iu -c. • - /!-_ . 


M * tf-j 1 ; 

v' 7 tr I -vi 

J . ) 


J 


J 
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The "dot" denotes different lation wLth respect to 'V'. It should 
be noted that INT^ - 0 and it can be shown that INT^ “ 

where ^ is the dissipation integral. The integrals INT^ and INT^ are 
moments of the Integrand of the dissipation integral. 

These are the equations programmed and solved in the present study. 

Once these equations are solved for the C.j, all of the boundary layer 
parameters of interest can be obtained by algebraic relations or staple 
quadratures. 

DISCUSSION OF METHOD 
General 

two general conclusions can be deduced from equation (7) . First, in 
application to the present formulation, it is clear that an arbitrary 
number of computing equations can be generated by choosing a suffta it number 
of linearly independent weighting functions. Therefore, solutions or an 
arbitrary number of coefficients in equation (8) m?\ be obtained. Second, and 
more generally, it can be seen that a whole family of computational methods 
are embodied in equation (7) . 'f for example, a large number of square 
wave weighting functions were employed, and some variation jf 8 with y were 
assumed over each increment, a finite difference formulation could be obtained. 
This correspondence between the generalized moment methods and finite difference 
, methods was recognized by Dorodnltsyn (1960) and also pointed out by Kendall 
and Bartlett (1967) . Similarly if a combined law of the wall and law of the 
wake were assumed for the velocity profile and f^(d) • 1, for 0 < ft < l, and 
fj(6) • 1 for 0 0 _< 1/2, a computing scheme similar to that employed by 
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Hoses (19WJ mid be obtained. 

In the present Method the coordinate vectors are chosen to be continuous 
over the entire region 0 < u ' 1 which alia! nates the need of establishing 
"internal* boundary conditions. 

The resulting systen of non-linear ordinary differential equations nay 
- ho integrated in the str ewr lse direction once the initial values of the C_% 
are chosen. This choice of the initial has proven to he a particularly 
difficult problem. To detezalne the Cj% representing experimental velocity 
profile at the initial station, too procedures mere considered. The first 
procedure employed was to Match the experimental values C^, {*, i, and kinetic 
energy thickness. <^. It uaa found that snail perturbations in the values 
of integral thickness parameters produced large eccursions in the values of 
the Gj. This difficulty led in a second proce du r e in which was matched and 
die predicted profile was made to pass through three data points in the 
initial experimental velocity profile. This procedure was found to yield 
satisfactory results for the integral parameters and was eaployed in the 
pre se nt study. Sines no clear rules for choosing the three velocity profile 
points have been established a trial and error procedure Bust be employed. 

To obtain Initial cc'dltlooa, it was found however that. Independent of the 
choice of points to be matched, any ant of yielding a good fit to the 
initial profile converged rapidly to a solution yielding die sane physical 
parameters. 

Choice of ) sad f^ 

The forms chosen for the weighting functions and the velocity gradient 
profile Oh are f'ven la the tceding section v Both of these choices were 
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w<l> on the basis of Che soccest. obtained using similar functions In 
lmdair flam analysis. The sue teas of the method depends on how well the 
approximate solution assumed represents the actual solution sought (i.e. , how 
rapidly the Infinite series representing the actual solution converges) and 
in aridities on how well the a ssu m e d solution can be resolved in the subspace 
- defined by the weighting fu nctio ns. This latter point has not previously 
been clearly enunciated by Investigators applying the method to boundary 
layer flows. In fact wher the function chosen for ) contains po wers of u not 
contained in the set of weighting functions, as In the present study sad in 
Pe i we rt and Abbott (1967) , an additional approximation is introduced. That is, 
the sesuaed solution does not lie wholly within the subspsce defined by the 
weighting functions. Since the solution to the ordinary differential equations 
is, in fact, obtained in tills subspace, only that portion of the initial 
profile which lies in the subspace nay be resolved in subsequent computations . 
When the representation for jl is so chosen, it is no longer possible to 
rigorously denons trace exactness and convergence of the general moment method. 

la tie present study the weighting functions contain terms up to fi J . 

It was found, h owe v er, that to represent the spectrum of profiles required in 
die present study, ana figure 1, sequentially complete poly n o mial s of order 
three were inadequate. In particular. It was found that In order to generate 
the sharply Inflected velocity profiles character is tic of favorable and zero 
pressure gradient flows large exponents must be chosen for the If representation 
(P » 2, Q ■ 4, E ■ 6). This choice of exponent was node for all flows 
represented in table 1 with the exception of the TUlnaim flow ID 1500, 
and the Hewn an flow ID 3500. For the Til Inarm flow. In order to obtain 
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correct Initial conditions , the exponent u h. ^ to be chosen as P - 1, Q ■ 2, 
1*3. In order to determine the effects of the choice of exponents in the 
) representation a parametric study was carried out using the Neuman flow 
as reference. Figure 2a shows a comparison of the results obtained for the 
exponents, ? ■ 1, Q • 5, 1*9 and the exponents ? » 1, Q - 2, R-3 with 
. die data of Hetman (ID 3500) . Five additional computer runs were made 
covering die Intervening range of exponents and it was found that, as the 
highest exponent decreased, at the last station, the skin friction coefficient 
decreased nonotonl cal 1,y , and the value of H Increased bopo topically . The 
nomen tun thickness, on die other hand, was essentially invariant with 
exponent chcic j and was in reasonable agreement with the data for all seven 
sets of computations. _It should be noted that die initial values (x^ » 2.759 
feet) of the integral parameters and skin friction were well represented by 
all die sets of exponents considered. It is clear from consideration of figure 
Zb that for the experimental velocity profile obtained by Neuman at station 
V, low order exponents most be employed. It is equally clear that none of 
the exponents considered yield profiles which adequately describe the exp«^i- 
nental results. Farther discussion of this point is deferred to the section 
entitled Discussion of Results. 

Choice of Eddy Viscosity Model 

The eddy viscosity nodel employed in this study is that of van Driest 
(1956) for the wall region and that of Cleuser (1956) for the wake region. 
These two models are simply patched together at equal values of c. A 
slight variation from tha van Driest model employed by Snlth et al. (1965) , 



Is used here in that the assuap i Ion of t * in chit sublayer has not been 
made in this stwiy. This model nay be written ss: 


» ^ *. t 


^ C.VC.Y* 


VJ 


c ■ 'fYzc.M \ \ 


The dauser Law Is simply' 

1 " — - V.C> -v O.C\C,S ^ 

M c 

This _aoice of eddy viscosity oodel was made simply to permit c li^iriaate 
comparison of the results of the present method with the exact solutions ox 
Smith ec al. (13i5). Such a comparison is shown in figure 3 for the flow 
of Schubaner and Klebanoff (ZD 2100). It can be seen that both calculation 
procedures yield essentially the same results. 

The basic moment method imposes no res trie tui ons on the choice of eddy 
viscosity nor does It require the explicit formulation of an eddy viscosity 
model. . This latter point is discussed in the following section. 


Numerical Solutions 

In an effort to make the present method comparable in machine time with 
the simpler integral methods an attempt was made to reduce the computation 
time involved in the streaavlae integration. * This can be accomplished by 
either reducing the number of integration steps required to traverse some fixed 
distance or by reducing the time required to take each step. Two methods of 
integration for the numerical solution of equation (10) were investigated and 
the time consumed by ev&Hacion if the dissipation Integrals was examined. 

The first integration procedure -mployed was fourth order variable 
step size Adams -Kaulcon predictor-corrector integration routine. Secondly. 
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a modified-Euler Implicit method employing local linearization of the 
derivatives was used. Lomax (1967) shows that the stability requirement 
Implicit in fourth order Adams-Moulton integration is [lb( < 0.7, where h 
is the Integration stepsize and A the largest eigenvalue of the locally 
linearised system. In order to avoid a small integration stepsize when 

parasitic eigenvalues are present (c.f. Lomax and Baily (1967)), the 
implicit Euler integration routine was employed. This latter integration 
scheme is stable for all | Ah| but requires local linearization of the 

3 

differential equations. The accuracy of the method is of 0(h) and the 

2 

error in the linearization is of course 0(h ) sc that some limitation of 
stepsize is required to retain sufficient accuracy. In order to control 
the error, the simple expedient of permitting a maximum 12 variation in the 
moat rapidly varying was employed. It was found that a factor of two 
reduction in machine time could easily be obtained using the implicit 
Integration. Table 1 compares IBM 7094 computation time for several sets 
of data for the two integration schemes. 

The second item considered was the reduction of time required for each 
step. It was found that roughly 70 percent of the time required for each 
step was associated with the evaluation of the dissipation integral and Its 
moments, (i.e. , the IMT^). Delwert and Abbott (1967) employed correlations 
of these integrals deduced from selected experimental data; Ttuckehbrodt 
(1955) , Goldberg (1966) , Rotta (1967) , and Nash and Spalding (given by Rotta) 
have each proposed "universal** correlations for the first of these Integrals. 
However, it is to be noted that these correlations differ substantially In 
both value and choice of correlation parameters. To employ this procedure 
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In the present four parameter furuulation requires correlations of tke 
dissipation Integral and its first and second moments. Since correlations 
of the first and second moments are not available in the literature, an 
effort was cade to determine approximate correlations of these quantities, 
the values of the three Integrals were computed from the eddy viscosity cod el 
for all flow conditions. Figure 4 shows the correlation of these results 
together with the correlation equations for the first integral given by 
Truckenbrodt (1955) end by Escudier and Spalding (1965) for the flat plate 
data of Vleghardt and TUlnaxm (ID 1400) . 

Substitution of these correlation equations into the computer prograc 
in place Of the nunerical evaluation of the integrals indicated chat the 
speed of computation could be increased by a factor of three but that the 
solutions appeared to he extremely sensitive to snail errors in the 
correlations . 

As s result of the above studies is Is concluded that if adequate 
correlations of the dissipation integrals were available and the implicit 

s 

Integration scheme were e — loyed a reduction of computer time of approximately 
an order of magnitude below those presented in table 1 for explicit integra- 
tion could be obtained with the present method. 

DISCUSSION OF RESULTS 
General 

The primary results presented in this paper are comparisons of predicted 
end experimental values of skin friction coefficient, end integral thickness 
parameters for sixteen different turbulent boundary layer flows. The 
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experimental flows considered act*, listed In cable 1 by identification ouster 
and the name(s) of the Investigdtor(s). More detailed results for representa- 
tive flow conditions together with an assessment of the adequacy of the 
present method are considered below. 

Experimental Data 

An examination of the data studied herein suggests a division into two 
distinct categories. This division is based on satisfaction of the streaevise 
integrated , two-dimensional momentum integral equation which may be written as: 

f d? v | 

J <** ^*e (11) 

X % ° 

Values of 6 and H (determined by integrating the experimental velocity profiles) 
and (determined by matching experimental velocity profiles and lav of the 
wall profiles) obtained at successive survey stations may be introduced 
into (11), along with the experimental valv s of C_ and dU^/dx, and integrated 
from x q to any downstream station x. Just such an analysis was carried out 
by Kline et al (1968) for all of the sets of data studied herein. It was 
found that for some flows, the left and right hand sides of equation (11) 
diverged at some streamwise station. Based on the above discussion, the data 
were divided into the two categories presented in table 2. The first category 
contains all of those flows, or portions of flows, which reasonably satisfy 
equation (11) and the second category contains those flows, or portions of 
flows, which do not satisfy equation (11). A detailed discussion of flows in 
the second category will not he given since any flows satisfying equations (1) 
and (2) must necessarily satisfy (11). However, comparisons of those data with 
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che present method will be presented for completeness. 

Comparison of Predicted and Experimental Results 

The comparisons of predicted and experimental values of H, and 
Reg are given in the figures showing the collected comparisons. For those 
flows satisfying the momentum integral balance, the predictions of B, C-, and 
Re^ are in good agreement with the data except for the flows with a large 
wake component. 

Figure 5 shows a comparison of experimental and predicted velocity profiles 
at four stations of the Schubauer and Klebanoff flow (ID 2100). As can be 
seen, good agreement between experiment and prediction is obtained everywhere 
accept at the 22 foot station where the boundary layer has been subjected to 
an adverse pressure gradient over the preceding five feet. 

Similar agreement at the initial stations and subsequent degradation 
of accuracy is shown in figures 6, 7, and 8 for the flows of Ludwieg and 
Tillmann (ID 1100), Claus er (ID 2200), and Ciauser (ID 2300), r - ^pectively. 

This lack of agreement gave rise to an investigation into the character of the 
velocity profiles which could be described by the currently employed velocity 
gradient profile, $. 

In an effort to determine the flexibility of the presently employed 
profile, a flow representative of these having large wake components, i.e., 
Ciauser (ID 2300) a^as examined. The velocity profile at the initial station 
was replotted on law of the wall coordinates in figure 9. It was found that 
the representation of the initial profile could be forced to agree with the 
data in the wake region, figure 9a, or in the wall region figure 9b, but not 
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In both regions. Therefore, It is concluded that the ^polynomial 
representation of the velocity profile is inadequate to represent velocity 
profiles with large wake cotaponents. 

In contrast to the above, flows with moderate wake components, e.g. Clauser 
(ID 2200), can be described reasonably well with the present formulation. 

■ This conclusion is substantiated by the results presented in figure 10. 

For those flows where the momentum integral balance is not satisfied, 
large wake components are usually present. This clouds the comparisons of 
predictions with data, since at the outset, equations (1) and (2) are invalid. 

COXCLUSIOKS 

The intent of the present study was to determine the efficiency and 
accuracy of a generalized moment method in the solution of the incompressible 
turbulent boundary layer equations. Relatively large computing times were 
required to obtain solutions for the flows considered. This was the result of 
the time consuming numerical evaluation of the dissipation integral and its 
two additional moments. By using an accurate correlation of these integrals 
the quoted run times can be reduced by approximately 70Z. The accuracy of 
the present formulation of the generalized moment method is acceptable for 
flows without large wake components. For flows with large wake components, 
the accuracy may be improved by employing more physically realistic velocity 
gradient profiles. 

It is believed that despite certain shortcomings of the present formula-* 
tlon, the moment method provides a powerful tool in the solution of partial 
differential equations. In addition, with a modification of the velocity 
profile to account for the presence of large wakes and the reduction in 
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computing time which has been sih.>wn to be possible, the net hod should -rove 
useful both in turbulent boundary layer research and as a design tool. 
Extension of moment methods, in particular, and any computing scheme, in 
general, to separating flows will require a substantially better understanding 
of the physical processe. involved. 
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TABLE I 


ID HUMBER 

INVESTIGATOR 

RES TIMES (7094 MIX.) 


- 

Explicit 

Implicit 

1100 

Ludvieg and Til Inarm 

4.77 


1200 

Loduieg and Tillaaim 

2.96 


1300 

Ludvieg and Tilluann 

5.03 


1400 

Ueghudt and Tilluann 

4.37 

1.85 

1500 

Tilluann 

10* 

2.95 

2100-I 1 

Schubaner am. KLebanoff 

10.97 

4.95 : 

2300-H 2 

Sebubauer tnd Ildaaoff 

2.00 


2200 

Clauser 

3.61 


. 2300 

Clauser 

5.94 


: 2400 

Bradshaw 

4.38 


2500 

Bradshaw 

1.61 


2600 

Bradshaw 

7.09 


2700-I 3 

Herring end Sorbury 

2.01 


2700-n 4 

Herring and Sorbury 

2.19 


2900 

Perry 

8.00 

3.70 

3500 

i Hewnan 

2.71 


3800 

Moses 

1.22 


4000 

Moses 

10* 

4.00 


1. » o - 1.5 feet 

2. x - 19 feet 
- o 

3. Investigators first pressure distribution. 

4. Investigators second pressure distribution. 
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(b) Matching velocity profile in vail region 

Figure 9 .- Comparison of initial velocity profiles for the present method 
with experimental data from a flow having a large wake component data 
of Clauser (ID 2300), x « 9.0 feet. 
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